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Abstract. The paper discusses the solution of a simple kinetic equation of 
the type used for the computation of the change of the chemical composition 
in stars like the Sun. Starting from the standard form of the kinetic equation 
it is generalized to a fractional kinetic equation and its solutions in terms 
of H-functions are obtained. The role of thermonuclear functions, which 
are also represented in terms of G- and H-functions, in such a fractional 
kinetic equation is emphasized. Results contained in this paper are related to 
recent investigations of possible astrophysical solutions of the solar neutrino 
problem. 
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1. Introduction 



A spherically symmetric, non-rotating, non-magnetic, self-gravitating model 
of a star like the Sun is assumed to be in thermal equilibrium and hydrostatic 
equilibrium, with a non- uniform chemical composition throughout. The star 
is characterized by its mass, luminosity, effective surface temperature, ra- 
dius, central density, and central temperature. For a given mass, four of 
these variables are independent and are governed by four simultaneous, non- 
linear, ordinary differential equations of the first order and four boundary 
conditions. Since there are four equations but more than four unknowns, 
additional information must be provided through the equation of state, nu- 
clear energy generation rate, and the opacity (constitutive equations). The 
assumptions of thermal equilibrium and hydrostatic equilibrium imply that 
there is no time dependence in the equations describing the internal structure 
of the star (Kourganoff, 1973, Perdang, 1976, Clayton, 1983). 

The evolution of a star like the Sun is governed by a second system of 
differential equations, the kinetic equations, describing the rate of change of 
chemical composition of the star for each species in terms of the reaction rates 
for destruction and production of that species (Kourganoff, 1973, Perdang, 
1976, Clayton, 1983). 

Methods for modeling processes of destruction and production have been 
developed for bio-chemical reactions and their unstable equilibrium states 
(Murray, 1989) and for chemical reaction networks with unstable states, oscil- 
lations, and hysteresis (Nicolis and Prigogine, 1977). Stabihty investigations 
of thermonuclear reactions of stellar interest have not yet been worked out in 
detail. However, the potentiality of instabilities in thermonuclear chains may 
not be overlooked, since, as was pointed out once by Eddington, "what is 
possible in the (Cavendish) Laboratory may not be too difficult in the Sun" 
(Perdang, 1976, Mestel, 1999). 

Consider an arbitrary reaction characterized by a time dependent quan- 
tity = N{t). It is possible to equate the rate of change dN/dt to a 
balance between the destruction rate d and the production rate p of N, that 
is dN/dt — —d + 'p. In general, through feedback or other interaction mecha- 
nisms, destruction and production depend on the quantity N itself: d — d{N) 
or p = p{N). This dependence is complicated since the destruction or pro- 
duction at time t depends not only on N{t) but also on the past history 



2 



N{t),t < t, of the variable N. This may be formally represented by 

dN/dt^-d{Nt)+p{Nt), (1) 

where Nt denotes the function defined by Nt{t*) = N{t — t*), t* > 0. Here d 
and p are functionals and eq. (1) represents a functional-differential equation. 
In the following we study a special case of this equation, namely the equation 

dN/dt = -aN{t) (2) 

with a constant a > 0. Eq. (2) imphes that spatial fluctuations or inho- 

mogenieties in the quantity N{t) are neglected. The standard solution of the 
differential equation (2) will be briefly discussed in Section 2 and the gener- 
alization to a fractional differential equation and its solution will be derived 
in Section 3. Conclusions will be drawn in Section 4. 



2. Standard Kinetic Equation 

The production and destruction of species is described by kinetic equations 
governing the change of the number density Ni of species i over time, that 
is, 

dN- 

-17^ = - E ^i^j < >ij + E ^kNi < av >ki, (3) 

j k,l^i 

where < av >mn denotes the reaction probability for an interaction involving 
species m and n, and the summation is taken over all reactions which either 
produce or destroy the species i (Haubold and Mathai, 1995). For a gas of 
mass density p, the number density iVj of the species i is expressed in terms of 
its abundance Xj, by the relation = pN^Xi/Ai, where Nj^ is Avogadro's 
number and Ai is the mass of i in mass units. The mean lifetime Tj{i) of 
species i for destruction by species j is given by the relation 

1 X ■ 

^^■(0 = —7^ = < av >ij^ pNa^ < TV >ij, (4) 

where is the decay rate of i for interactions with j. Eq. (4) reveals the 
physical importance of < av >ij for the kinetic equation (3). 
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In the case of a nondegenerate, nonrelativistic gas, if a nonresonant 
charged nuclear reaction proceeds at low energies dominated by Coulomb- 
barrier penetration, the reaction probability < av >mn takes the form (Clay- 
ton, 1983, Bergstroem et al., 1999) 



(5) 



where I2 represents a thermonuclear function given by 



hiiy -l,a,z,p) = / dyy'"^e-°-y-^y ' , i/>0,a>0,^>0,p>0 

Jo 



a 0-2,0 



1 




azp 


(^,i),(o,^) 



(6) 



where H denotes Fox's H-function (Mathai, 1993, Haubold and Mathai, 
1998), which was introduced into mathematical analysis to unify and extend 
existing results on symmetrical Fourier kernels. 

Fox's H-function (Mathai and Saxena, 1973, 1978) is defined in terms of 
a Mellin-Barnes type integral 



where 



H{z) 
h{s) = 



Tjm.n 



|(ai,ai),...,(ap,ap) 
l(6i,/3i),...,(6<,A) . 



1 



h{s)z ^ds, 



2Tii Jl 

{n- r(/;, + /i,,s)}{n" r(i-«,-n,,s)} 



(7) 
(8) 



i — V— 1, and L is a suitable path which will be briefly described in the follow- 
ing. An empty product is interpreted as unity and it is assumed that the poles 
oir[bj+f3js), j = l,...,m, are separated from the poles of r(l— a^— a^s), j = 
1, . . . ,n;ai, . . . ,ap,bi, . . . ,bq are complex numbers; ai, . . . , ap, . . . , /3g are 
positive real numbers. The poles of T{bj + /3js),j = 1, . . . ,m, are at the 
points 

s = -(bj + v)l(5j, j = 1, . . . , m, i/ = 0, 1, . . . 
and the poles of r(l — aj — CKjs),j — 1, . . . ,n, are at 

s (1 - Ofc + \)/ak, = 1,2, . . . ,n, A = 0, 1, . . . 

The condition of separability of these two sets of poles imposes that there be 
a strip in the complex s-plane where the H-function has no poles. There are 
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three types of paths L possible. These correspond to paths 1,2, and 3 shown 
in Fig. 1. 



Fig. 1. 

For all practical problems where H-functions are to be applied one mainly 
requires paths 2 and 3. Hence the conditions coming from these two paths 
are given in the following. It is to be pointed out that when more than one 
path L makes sense then it can be shown that they lead to the same function 
and thus there will be no ambiguity. 
Let 

/^ = E/5.-E«. (9) 

and 

P = Ii%,afTlUp/\ (10) 
The H-f unction exists for the following cases: 



Case (i) g > 1, > 0, H{z) exists for all z^z ^ 

Case (ii) q> = 0,H(z) exists for \z\ < f3~^. 

Case (iii) p > 1, < 0, H{z) exists for all z,z 0. 

Case (iv) p > 1, // = 0, H{z) exists for \z\ > P~^. 



In the above cases it is assumed that the basic condition is satisfied that 
the poles of T(bj + Pjs),j = 1, . . . ,m and r(l — aj — ajs),j = 1, . . . ,n are 
separated. Note that in cases (i) and (ii) the H-function is evaluated as the 
sum of the residues at the poles of T(bj + Pjs),j = 1, . . . , m and in cases (iii) 
and (iv) the H-function is evaluated as the sum of the residues at the poles 



of rri 



When ai 



ajs),j 



.n. 



Pi — . . . — Pq = 1, the H-function reduces to 



Meijers's G-function. That is. 



l(6i,l),. 



Tjm.n 



■ ,{ap,l) 



P,<1 



\bi, 



(11) 



When p in cq. (6) is real and rational, then the H-function can be reduced 
to a G-function by using the multiplication formula for gamma functions. 



r(mz) = (27rY-^m'^'-^U"K}r ( z + ^] ,m ^ 1,2, . . . 

■' \ mj 



(12) 
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In the case of Coulomb-barrier penetration (Gamow factor), p = \il2 reduces 
to 



T ( ^\ " Vv3,0 
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(13) 



where G denotes Meijer's G-function, which was introduced into mathemat- 
ical analysis in attempts to give meaning to the generalized hypergeometric 
function pFq in the case p> g + 1. 

Proceeding with eq. (3), the first sum in eq. (3) can also be written as 

- NiNj < av >ij^ -Ni(J2 < >ij) = (14) 

j j 

where Oj is the statistically expected number of reactions per unit volume 
per unit time destroying the species i. It is also a measure of the speed in 
which the reaction proceeds. In the following we are assuming that there 
are Nj{j — 1, . . . ,i, . . .) species j per unit volume and that for a fixed Ni 
the number of other reacting species that interact with the i-th species is 
constant in a unit volume. Following the same argument we have for the 
second sum in eq. (3) accordingly, 

+ ^ NkNi < av >ki^ +Nibi, (15) 

where Nibi is the statistically expected number of the i-th species produced 
per unit volume per unit time for a fixed Ni. Note that the number density 
of species i, = Ni{t), is a function of time while the < av >mm containing 
the thermonuclear functions (see eqs. (5) and (6)), are assumed to depend 
only on the temperature of the gas but not on the time t and number densities 
Ni. Then eq. (1) imphes that 

^ = -(a,-W(t). (16) 

For eq. (16) we have three distinct cases, Ci = ai — bi > 0, q < 0, and q = 0, 
of which the last case says that Ni does not vary over time, which means 
that the forward and reverse reactions involving species i are in equilibrium; 
such a value for Ni is called a fixed point and corresponds to a steady-state 
behavior. The first two cases exhibit that either the destruction (q > 0) of 
species i or production (q < 0) of species i dominates. 
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For the case q > we have 

dt 



= -CiNiit), (17) 



with the initial condition that Ni(t = 0) — Nq is the number density of 
species i at time t — 0, and it foUows that 

Ni{t)dt = Noe-^'^dt. (18) 

The exponential function in eq. (18) represents the solution of the linear 
one-dimensional differential equation (17) in which the rate of destruction of 
the variable is proportional to the value of the variable. Eq. (17) does not 
exhibit instabilities, oscillations, or chaotic dynamics, in striking contrast to 
its cousin, the logistic finite-difference equation (Perdang, 1976, Haubold and 
Mathai, 1995). A thorough discussion of eq. (17) and its standard solution 
in eq. (18) is given in Kourganoff (1973). 



3. Fractional Kinetic Equation 

In the following, for the sake of brevity, the index i in eq. (17) will be 
dropped. The standard kinetic equation (17) can be integrated 

N{t)-No = -CoD;^N{t), (19) 

where qD^^ is the standard Riemann integral operator. The generalization 
of this operator to the fractional integral of order p > is denoted by aDt^ 
and is defined, following Riemann- Liouville, based on the Cauchy formula, 
by 

aA"7(t) = [ drf{T){t - Tf-\ p > (20) 

with 

aD'jit) = fit) 



(Oldham and Spanier, 1974, Miller and Ross, 1993). The most general frac- 
tional integral operator of the type (20) contains Fox's H-function (7) as the 
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kernel function. If f{t) is continuous for t > a, then integration of arbitrary 
real order has the property 



Replacing the Riemann integral operator by the fractional Riemann-Liouville 
operator qD^" in eq. (19), we obtain a fractional integral equation corre- 
sponding to eq. (17) 



(21) 



For dimensional reasons, the coefficient c in eq. (19), containing the probabil- 
ities of the reaction under consideration, had to be replaced by c'^ accordingly. 
The Laplace transform of the Riemann-Liouville fractional integral is 



where 



L{oD^Pf{ty,p}=p-PF{p), 
F(p) = r(p) / dre-P-fir). 

Jt=0 



(22) 



In order to solve eq. (21), the integral equation is exposed to a Laplace 
transformation leading to 



N{p)^L{N{ty,p} = No 



P' 



l + (f)-'^' 



(23) 



To arrive at a representation of eq. (23) in terms of Fox's H-function, 
Mathai and Saxena's (1978) result can be used. 



y0 



1 + az"" 
N{p) 



a « III 1 



No-Hl:l 
c 





V 

I 











(24) 



To prepare eq.(24) for an inverse Laplace transform, the following two fun- 
damental properties of an H-function can be used. 



)^ P'<1 



,(ai,ai),...,(ap,ap) 
l(bl,/3i),...,(b„/3,) 



k I {ai,kai),...,{ap,kap) 
" l(6i,fc/3i),...,(6„fe/3,) 



, A; > 0, (25) 



Tjm.n 



<{ai,ai),...,(ap,ap) 
l(&i,/3i),-,(6<,A) 



Tjn.m 



l(l-ai,ai),...,(l-ap,ap) 



leading to 



1 



N{p) = No—Hli 



-p 






.c 







(26) 
(27) 



where the H-function is defined in eq. (7). 

The Laplace transform of Fox's H-function (7) is given in terms of another 
H-function by 



L{H{zy,p}^^H^;il'r 



p 



(l,l),(l-ap,ap) 



(28) 



for < < 1 in (9), and 



rn,n+l 



L{H{zy,p} = -H;rc 



(0,l),(ap,ap) 



(29) 



for /X > 1 in (9), respectively. 

Further, having H{p), the inverse Laplace transform of this H-function is 
given by 

1 

z 

z 



H{z) = L-'{H{p),z} = -H-;;i, 
for < // < 1 in (9) and 



(i-b.A) 

(l-ap,ap),(l,l) 



(30) 



H{z)^L-'{H{p),z}^-H;:r,,^ 



(ap,ap), (0,1) 



(31) 



for > 1 in (9), respectively. 

The above four Laplace transforms hold for 

maxi<j<„ 3? < mini<j<„ 3? . 

Applying an inverse Laplace transform to the H-function in eq. (27) gives 



iV(t) = N^-E]:l 



ct 



(0,-i) 

{0,|),{0,1) 



(32) 



which is the solution of the fractional kinetic equation (21). For the H- 
function in eq. (7) with (32), the following computable representation can 
be derived (Mathai and Saxena, 1978). When the poles of Wj^^-^r{hj — (3js) 
are simple, that is, 

I3h{hj + \)^l3,{hh + v) 
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for j 7^ h;j,h = 1, . . . ,nr, X, u = 0, 1, 2, ... . Then one obtains the following 
expansion for the H-function, 



tlto {^=m+l^ {l - + 

{n^..,,r{a,-a,i^}} ' ^ ^ 

which exists for all 2; 7^ if /x > and for < \z\ < if = 0, where 
^ and (3 are given in eqs. (9) and (10). Comparing (32) with eq. (33), one 
obtains the series expansion 

^W = /V„gf(Ll_M-. (34) 

For ly — 1, the exponential solution of the standard kinetic equation (18) is 
recovered. 



4. Conclusions 

The gravitationally stabilized solar fusion reactor, because of its density and 
temperature, is not a weakly interacting or high temperature plasma. Re- 
cently, Kaniadakis et al. (1997) and Coraddu et al. (1998) have explored the 
possibility that the electrical microfield distribution influences the particle 
dynamics, collisions are not two-body processes and retain memory beyond 
single scattering events, and that the velocity correlation function has long- 
time memory arising from the coupling of collective and individual degrees 
of freedom. In this connection particle diffusion is a non-Markovian pro- 
cess (anomalous diffusion) and diffusion and frictional coefficients are energy 
dependent. Based on these results, they conclude that the equilibrium sta- 
tistical distribution function differs from the Maxwell-Boltzmannian one and 
is governed by generalized Boltzmann-Gibbs statistics. 

Even if the deviations from the Maxwell-Boltzmann distribution, that are 
compatible with the current knowledge of the solar core behavior, are small, 
they arc sufficient to sensibly modify the sub-barrier particle reaction rates 
and subsequently solar neutrino fluxes. The above authors have also shown 
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that the respective modifications of the reaction rates do not affect bulk 
properties of the gravitationally stabihzed solar core such as sound of speed 
or hydrostatic equilibrium which depend on the mean values obtained by 
averaging over the Maxwell- Boltzmann distribution function (Degli'Innocenti 
et al. 1998). 

Haubold and Mathai (1998) have derived general closed-form represen- 
tations of particle reaction rates that are suitable to incorporate changes of 
the Maxwell-Boltzmann distribution function, mindfully taking into account 
the fact that the whole distribution has physical meaning, contrary to the 
case of, say, the Gaussian law of errors. In this paper we proceeded one 
step further in generalizing the standard kinetic equation (17) to a fractional 
kinetic equation (21) and derived solutions of a fractional kinetic equation 
that contains the particle reaction rate (or thermonuclear function) (5) as 
a time constant, and provided the analytic technique to further investigate 
possible modifications of the reaction rate through a kinetic equation. The 
Riemann-Liouville operator in the fractional kinetic equation introduces a 
convolution integral with slowly-decaying power-law kernel, which is typical 
for memory effects referred to in Kaniadakis et al. (1997) and Coraddu et 
al. (1998). This technique may open an avenue to accommodate changes in 
standard solar model core physics (Dar and Shaviv, 1999) as proposed by 
Shaviv and Shaviv (1999) and Schatzman (1999). 

In the solution of the fractional kinetic equation (21), given in eqs. (32) 
and (34), the standard exponential decay is recovered ior u — 1. However, 
eqs. (32) and (34), for < z/ < 1, show power law behavior for i — > oo and 
are constant (initial value Nq) for t ^ 0. In the investigations in this paper, 
we used, as an example, the standard thermonuclear function 72°°'' = Ii, the 
same computations can be done for l2'^\ J3, or J4 with additional parame- 
ters showing modifications of the Maxwell-Boltzmann distribution function 
(Haubold and Mathai, 1998). 

All analytic results in this paper have been achieved by the application 
of the theory of generalized hypergeometric functions, particularly Meijer's 
G-function and Fox's H-function, which seem to be the natural means for 
tackling the problems referred to above. 
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